/*
 * Reaction.cpp
 *
 *  Created on: 21 Jul 2011
 *      Author: Allan
 */

#include "Reaction.h"

// C++ includes
#include <cmath>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "NumericalUtilities.h"

Reaction::Reaction()
{}

Reaction::~Reaction()
{}

void Reaction::SetReactant(const Species& reactant)
{
	this->reactant = reactant;
}

void Reaction::SetProducts(const vector<Species>& products)
{
	this->products = products;
}

void Reaction::SetEquilibriumConstantData(const EquilibriumConstantData& eqData)
{
	this->eqData = eqData;
}

void Reaction::SetSpeciesIndexes(const SpeciesMap& speciesMap)
{
	// Set the index of the reactant
	reactant.index = speciesMap[reactant.name];
	
	// Loop over the products of the reaction and set their indexes
	BOOST_FOREACH(Species& product, products)
		product.index = speciesMap[product.name];
}

const Reaction::Species& Reaction::GetReactant() const
{
	return reactant;
}

const vector<Reaction::Species>& Reaction::GetProducts() const
{
	return products;
}

const Reaction::EquilibriumConstantData& Reaction::GetEquilibriumConstantData() const
{
	return eqData;
}

const double Reaction::EquilibriumConstant(double T, double P) const
{
	// The base 10 logarithm of the equilibrium constant, i.e. log(K) 
	const double log10_K = LagrangeInterpolation(T, &eqData.Ti[0], &eqData.logKi[0], eqData.Ti.size());
	
	// Return K, not log(K)
	return std::pow(10.0, log10_K);
}

const double Reaction::ReactionQuotient(const DoubleVector& a) const
{
	// The reaction quotient of the reaction
	double Qr = 1.0/a[reactant.index];
	
	// Loop over all the products in the reaction
	for(uint i = 0; i < products.size(); ++i)
	{
		const double ai = a[products[i].index];      // The activity of the product i
		const double vi = products[i].stoichiometry; // The stoichiometry of the product i
		
		Qr *= std::pow(ai, vi);
	}
	
	return Qr;
}

/// --- Utility methods --- ///
const MatrixXd AssembleReactiveMatrix(const vector<ReactionPtr>& reactions, const StringVector& primarySpecies)
{
	// The number of reactions (Nr) and primary species (Nj)
	const uint Nr = reactions.size();
	const uint Nj = primarySpecies.size();
	
	// The stoichiometric matrix for the set of reactions
	MatrixXd vrj(Nr, Nj);
	
	// Loop over all the reactions
	for(uint r = 0; r < Nr; ++r)
	{
		// Loop over all the primary species
		for(uint j = 0; j < Nj; ++j)
		{
			// Get a reference to the vector of product species of the current reaction 
			const vector<Reaction::Species>& products = reactions[r]->GetProducts();
			
			// Find the (local) index of the product that corresponds to the current primary species j (note that the current reaction may not have this primary species as a product)
			const uint iproduct = std::find_if(products.begin(), products.end(),
					                  [&](const Reaction::Species& product) { return product.name == primarySpecies[j]; }) - products.begin();
			
			// Check if iproduct is within the limits [0, numProducts)
			vrj(r, j) = (iproduct < products.size()) ? products[iproduct].stoichiometry : 0.0;
		}
	}
	
	return vrj;
}

ostream& operator<<(ostream& out, const Reaction& reaction)
{
	out << reaction.GetReactant().name << " (" << reaction.GetReactant().index << ")" << " <=> ";
	
	BOOST_FOREACH(const Reaction::Species& product, reaction.GetProducts())
		out << setprecision(2) << product.stoichiometry << " " << product.name << " (" << product.index << ")" << ", ";
	
	return out;
}
